On steady-state currents through nano-devices: a scattering-states numerical 
renormalization group approach to open quantum systems 
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We propose a numerical renormalization group (NRG) approach to steady-state currents through 
nano-devices. A discretization of the scattering-states continuum ensures the correct boundary 
condition for an open quantum system. We introduce two degenerate Wilson chains for current 
carrying left and right-moving electrons reflecting time-reversal symmetry in the absence of a finite 
bias V. We employ the time-dependent NRG to evolve the known steady-state density operator 
for a non-interacting junction into the density operator of the fully interacting nano-device at finite 
bias. We calculate the differential conductance as function of V, T and the external magnetic field. 
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Introduction: The description of quantum systems 
out of equilibrium is one the fundamental challenges in 
theoretical physics. Even a simple non-equilibrium situ- 
ation, the current transport through an interacting junc- 
tion at finite bias is not fully understood. The Coulomb 
blockade [l[ and advent of the experimental realizations 
of the Kondo effect in such devices 0,0] requires a many- 
body description at low temperatures. 

While the equilibrium dynamics is well understood 
the non-equilibrium steady-state has been mainly investi- 
gated using perturbative approaches [M JH, [3 , IH based on 
Keldysh theory Q, the Toulouse point [l0| and the flow 
equation 11 1 . Landauer-Buettiker type approaches [l^ 
treat the charging effect only on a mean-field level by 
mapping the strongly interacting quantum problem onto 
a model of non-interacting fictitious particles, unsuitable 
to describe the Coulomb-blockade physicsfl]. In weak 
coupling and high temperature, the ac and dc-transport 
through molecular wires can be addressed by a quantum 
master equation for the reduced density matrix of the 
i unction [ 1 31] . All those approaches have only a limited 
validity of their parameter regimes. Recently, Han pro- 
posed an alternative perturbative metho d |1 4| based on 
Hershfield's steady-state density operator Tail, 12, HI- 
Based on similar ideas, a scattering-states Bethe-ansatz 
approach to an interacting spinless quantum dot has been 
implemented (l9| for finite bias. 

We present a numerical renormalization group 
approach [4| t o open quantum systems based on scatter- 
ing states [1 5]. It combines (i) Wilson chains for single- 
particle scattering states proposed below, (ii) Hershfield's 
steady-state density operator [3] for a non- interacting 
junctions at finite bias and (iii) the time-dependent NRG 
(TD-NRG)0, ED, 0]. Our scattering-states basis will 
be also useful for Quantum Monte Carlo and density 
matrix renormalization group (DMRG) approaches [23|. 
With our non-perturbative method, steady-state currents 
through interacting nano-devices can be obtained accu- 
rately for arbitrary temperatures, magnetic fields and in- 
teraction strength. 



Dissipative steady-state currents only occur in open 
quantum system in which the system size L has been 
sent to L — > oo before t — > oo. Transient currents can be 
calculated on a finite-size tight-binding chain within the 
TD-NRG as well as the time-dependent DMRGfH, Hi]. 
However, such transient currents vanish for t — > oo or 
even reverse their sign[24| in those approaches, a con- 
sequence of the non-interchangeable limit t — > oo and 
L — > oo[l8(. We circumvent this problem by discretizing 
a single-particle scattering states basis. Therefore, those 
states remain current carrying and a faithful representa- 
tion of an open quantum system. 

Theory: Interacting quantum dots (QD), molecular 
junctions or other nano-devices are modelled by the in- 
teracting region Hi mp , a set of non-interacting reser- 
voirs H.B and a coupling between both sub-system TLj: 
T~i = Ti-imp+Ti-B+'Hi ■ Throughout this paper, we restrict 
ourselves to a junction with a single spin-degenerate or- 
bital d with energy E^, subject to an external magnetic 
field H and an on-site Coulomb repulsion U. The or- 
bital is coupled to a left (L) and a right (R) lead via the 
tunneling matrix elements V a =L,Ri and T~C given by 
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Here h d a = d\d a , and c\ aa creates a conduction electron 
in the lead a of energy e and density of states p(e). 

This Hamiltonian is commonly used to model ultra- 
small quantum dots[H, In the absence of the local 
Coulomb repulsion Hu = U(J2 a nt - !) 2 / 2 , t hc 

single 

particle problem is diagonalized exactly in the contin- 
uum limit (l4 . 15 . T5 17 . 25|, [26| by the scattering states 
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FIG. 1: The local d-orbital is expanded in left-moving and 
right-moving scattering states. Each contributions defines one 
fictitious local orbital d aa of the junction of the scattering- 
states NRG. The Coulomb repulsion introduces backscatter- 
ing between left and right-movers. 
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where V — \/V£ + V R , and the Green function Gq ct (;z) = 

[z - (E d + U/2 - <tH/2) -V 2 J dep(e)/(z - e)] ~\ In 
the limit of infinitely large leads, the single-particle spec- 
trum remains unaltered, and these scattering states di- 
agonalize the Hamiltonian [l5| (p} for U = 0: 



Hl=H(U = 0) = ]T I d ^l a % 
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Hershfield has shown that the density operator for such a 
non-interacting current-carrying quantum system retains 
its Boltzmannian form (TBI fla ] even at finite bias: 



Po 



Tr 



The Yq operator accounts for the occupation of the left 
and right-moving scattering states, and p a for the differ- 
ent chemical potentials of the leads. 

Steady-state NRG: In order to apply the NRG to such 
an open quantum systems, the scattering states 7 ect(T are 
discretized on a logarithmic energy mesh using the NRG 
discretization parameter A0|. In contrary to a closed 
system, however, each of these single-particle states car- 
ries a finite current. Even for asymmetric coupling, the 
spectra of the right and left-movers remains symmetric, 
and the total current vanishes always at zero bias. 

Defining the creation operator for a fictitious left or 
right-moving c^Q-orbital d\ a — V J de^J p(e)[Gl a (e + 
^)]*7e<ra! the physical d- level can be decomposed into 
d\ = r R d) aR + Tl($ L by inverting Eq. |(3J) and using 
fa = V a /V. For U = 0, the Hamiltonian is diagonal 
in the left and right-movers. We use these d aQ -orbitals 
as starting vector of the Householder transformation [1| 
mapping the discretized scattering states continuum onto 
two semi-infinite Wilson chains[4|, as depicted in Fig. [T] 
These chains are almost identical to standard Wilson 



chain of a non-interacting resonant level model [4(. Each 
fictitious dg-Q-orbital consists of a normalized linear com- 
bination of scattering states 7 CCTQ ,: no auxiliary degrees 
of freedom has been introduced into the problem! 

We divide G r 0tJ (e + iS) into magnitude and phase, 
GQ a (e + iS) = e 1 ®"^ |Gg CT (e + iS)\, and absorb the energy 
dependent phase $<r(e) into the scattering-states opera- 
tors 7 C(T q by a gauge transformation. Then, the Wilson 
chains consist only of purely real tight-binding parame- 
ters. Diagonalizing the proposed scattering-states Wil- 
son chains yields a faithful representation of the steady- 
state density operator p for arbitrary bias. 

The current operator expanded in scattering states 
Itaa acquires an additional energy dependence via the 
scattering-phase shift $ CT (e). In our model (JTJ) , however, 
the current remains connected to the spectral function 
Ad(u>) of the retarded non-equilibrium Green function 27 1 



s~t roo 

I{V) = — E / duj M w " - /(« ~ Mil)] 

XTTA da (u)T (5) 

in such a scattering-states formulation even for finite 
U\ui E3, EH- /H denotes the Fermi function, G = 
(e 2 /h)4T L T R /(T L + r R ) 2 , T a = r 2 a nV 2 p(0), T = T L +T R 
and irA da {u}) = —^smG r da {u) + iS). 

Coulomb interaction: Expanding the operator n d a in 
the orbitals d aa yields two contributions: a density and 
a backscattering term: n d a = ffi a + 0„ , with n° a = 
12a r a4a^a- The backscattering term reads 
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and describes transitions between left and right-movers. 
This term vanishes in the tunnelling regime, where cither 
Tjj or tr vanishes. 

We will include the full Coulomb interaction into our 
theory in two steps. Since H°j, defined as = 
■j (^ CT — l) 2 , commutes with Yq, the steady state 

density operator po evolved into po — exp[— (3(H l — Yq)]/Z 
with Ji 1 = H + proven by the arguments given in 
Ref. [lH . 0^ ack can be neglected in the tunneling regime 
where p — > p . Then, the steady-state spectra is com- 
pletely determined by a single effective orbital, and the 
equilibrium spectral function is recovered. 

Ti. 1 marks the new starting point of our theory. The 
full Hamiltonian Ti. of the interacting model differs from 
TL 1 by the additional backscattering terms. TL does not 
commute with Yq, and the analytical form of steady-state 
density operator of the fully interacting problem is not 
explicitly known 15, 18}. We obtain a solutionis 21, 22j 
by evolving po with respect to the full Hamiltonian Ji into 
its steady-state value p^ = lim^oo e~ lUt p e lHt . In the 
current- voltage relation (0, the spectral function Ad a {<jj) 
for U = is replaced by the non-equilibrium spectral 
function[28| calculated with respect to p^. The details of 
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FIG. 2: (color online) Non-equilibrium spectral function for 
(a) a symmetric junction R = 1 at various values of finite bias 
voltage V, and (b) for a strongly asymmetric junction R = 
1000. The insets show the evolution of the Kondo-resonance. 
Parameters: U = 8, e/ = —4 and T — » 0. 
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this algorithm embedding the calculation of equilibrium 
spectral functions [29j, |3fJ are published in Ref . [22| . 

Results: All energies are measured in units of 
r = nV 2 p(0); a constant band width[4] of p(u>) = 
1/(2D)Q(D - \u)\) was used with D/T = 10. The num- 
ber of retained NRG states was N s = 2200; a A = 4 
was chosen. The model lacks channel conservation: only 
the total charge and z-component of the spin served as 
quantum numbers. We defined R = Tl/Tr and always 
kept r = Tl + Tr constant. The two chemical potentials 
fi a were set to fiL — —"^bY anc ^ ~ "^lY as function 
of the external source-drain voltage V consistent with a 
serial resistor model. 

The non-equilibrium spectral function for a symmetric 
junction is plotted for U — 8 and different bias V in 
Fig. [Ha). Multiple backscattering events cause gain (or 
lost) of single-particle excitation energy proportional to 
the applied bias. The Kondo resonance is destroyed with 
increasing bias due to redistribution of spectral weight 
towards higher energys. An onset of two weak peaks 
in the vicinity of the two chemical potentials remains 
for |V| > r[16|. For large R 3> 1 such backscattering 
processes are suppressed. The spectral function remains 
bias-independent. The Kondo resonance remains pinned 
to hl — > as depicted in Fig. [2](b) , and we recover the 
tunneling regime. 

The differential conductance is plotted for different 
asymmetry ratios R in Fig. [3Ja) using the same parame- 
ters as in Fig. O With increasing R, the non-equilibrium 



FIG. 3: (color online) The differential conductance G = 
dl/dV as function of the bias voltage (a) for different 
asymmetry factors R, (b) for different magnetic field H — 
0,0.1,0.2,0.4 and R = 1. Parameters: as in Fig. (c) Com- 
parison between the results for U = 5 from Ref. [l6l ] and the 
NRG calculation at T/F = 0.04 and R = 1 using z-averaging 
over 4 z-values[H 

spectral function is less broadened and, therefore, G(V) 
decreases for large bias voltage. Asymptotically, G ap- 
proaches the equilibrium t-matrix which is the exact re- 
sult for R — > oo and T —> 0. 

The effect of an external magnetic field onto the differ- 
ential conductance is shown in Fig. [3^b) . An increasing 
magnetic field splits the zero-bias anomaly which is fur- 
ther suppressed by the finite bias in a symmetric junction. 
This field dependence has been used in experiments [2| as 
hallmark for the Kondo physics at low temperatures. 

In Fig. (3jc), the NRG conductance for U = 5 is com- 
pared to the result of Ref. [l6| . Both curves agree for low 
bias. The NRG result shows a weaker decay of the zero- 
bias anomaly with increasing bias with a less pronounce 
maximum at large bias. The symmetrized equilibrium 
t-matrixQ is added for comparison as dashed line. 

The more generic case of an asymmetric junction with 
respect with a relatively large local Coulomb repulsion 
is plotted in Fig. |U The differential conduction reflects 
the lack of symmetry under source-drain voltage reversal. 
As depicted, the zero-bias peak vanishes with increasing 
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FIG. 4: (color online) The differential conductance G as func- 
tion of the bias voltage for different temperatures. Parameters 
R = 10, e/ = -1.5 and U = 12. 



temperature. 

Conclusion: A powerful new approach to the steady- 
state currents through nano-devices has been presented. 
We have introduced a NRG method based on scatter- 
ing states to incorporate the correct steady state bound- 
ary condition of current carrying systems. The steady- 
state density operator 15j of a non-interacting junction 
is evolved into the one of the interacting nano-device 
using the TD-NRG 20]. We have established an accu- 
rate solution for the strong-coupling regime and calcu- 
lated steady-state currents for arbitrary ratios R at fi- 
nite bias. The tunneling regime is included as an ex- 
act limit. Our approach does not suffer from any cur- 
rent reflection inherent to numerical simulations of closed 



quantum systems [24j. We have concentrated on the low- 
temperature properties of the nano-device, since the com- 
bination of arbitrary bias, large Coulomb repulsion and 
finite magnetic field remains the most difficult regime for 
all perturbative methods. However, the NRG is equally 
suitable to calculate the crossover from the low to the 
high-temperature regime as demonstrated in Fig. 01 An 
experimental hallmarkQ for Kondo physics, the splitting 
of the zero-bias Kondo peak with increasing magnetic 
field, is correctly described by our approach for arbitrary 
temperature, bias and field strength. 

This theory can be extended to more complicated 
multi-orbital models. Eq. (O must be modified and re- 
quires more complex correlation functions. Since single- 
particle scattering states can always be obtained ex- 
actly, the construction of the Wilson chain parameters 
is straight forward using the corresponding expansion of 
the local degrees of freedom and combining it with the 
transformation used for non-constant density of states 
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